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Abstract — Waveguide and resonant properties of diffractive 
structures are often explained through the complex poles of their 
scattering matrices. Numerical methods for calculating poles of 
the scattering matrix with applications in grating theory are 
discussed and analyzed. A new iterative method for computing 
the scattering matrix poles is proposed. The method takes account 
of the scattering matrix form in the pole vicinity and relies 
upon solving matrix equations with use of matrix decompositions. 
Using the same mathematical approach, we also describe a 
Cauchy-integral-based method that allows all the poles in a 
specified domain to be calculated. Calculation of the modes of 
a metal-dielectric diffraction grating shows that the iterative 
method proposed has the high rate of convergence and is 
numerically stable for large-dimension scattering matrices. An 
important advantage of the proposed method is that it usually 
converges to the nearest pole. 

Index Terms — quasiguided eigenmode, optical resonance, scat- 
tering matrix pole, diffraction grating 

I. Introduction 

DIFFRACTIVE micro- and nanostructures with resonant 
properties are of great interest when designing mod- 
ern elements of integrated guided-wave optics and photonics 
(photonic crystal fibers, waveguide grating couplers, optical 
sensors, guided-mode resonant filters, laser resonators) [[T)- 
|14|. Appearing as an abrupt change in the transmittance 
and reflectance spectra, the resonant and waveguide properties 
are normally associated with the excitation of the structure's 
eigenmodes. The structure's eigenmodes can be described 
in terms of the complex poles of the scattering matrix 0, 
0. Such an approach to explaining optical properties of 
the diffractive structures has been widely employed when 
describing the optical properties of diffraction gratings (2D 
and 3D) H, 0, including those comprising anisotropic lfT31 
and gyrotropic [9| materials, photonic crystal structures lT72l . 
Ifl3l and laser resonators lfl4l . 

The important practical problem requiring calculation of the 
poles of the scattering matrix is the design of grating/photonic- 
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crystal waveguides with prescribed optical properties (with 
specified dispersion of waveguide's quasiguided or leaky 
modes). This design problem belongs to the class of inverse 
problems and is solved with use of specialized optimization 
techniques. When solving this inverse problem, the compu- 
tationally effective methods for solving the forward problem 
(i.e. calculation of the scattering matrix poles) are of great 
practical importance. 

The calculation of the scattering matrix poles is computa- 
tionally challenging. A number of methods for solving this 
problem have been proposed in recent papers. The simplest 
techniques calculate the poles of the scattering matrix deter- 
minant lfl2l . Ifl3l or poles of its maximal eigenvalue |16|. 
A more advanced method proposed in Refs. Q, (8), ifTTl 
relies upon the linearization of the scattering matrix inverse. In 
order for the optical properties of metal-dielectric diffraction 
gratings and plasmonic structures to be adequately described, 
the scattering matrix dimension should be rather high. Usu- 
ally, the scattering matrix is calculated using Fourier modal 
method fl8l . Despite a number of approaches developed to 
enhance the method convergence 02), fl9l . the scattering 
matrix dimension can still amount to several hundreds [20 1 . 
For scattering matrices of such large sizes operations of 
calculating the determinant and the scattering matrix inverse 
in Refs. Q, 0, C3, D3 become numerically unstable ED, 
thus essentially limiting the methods' applicability area. 

Of special note is the method that calculates the eigenmode 
frequencies using the Cauchy integral. This method is used 
when calculating modes propagating in a slab waveguides [0, 
O and photonic crystal structures [6 |, [16], [22]. It is notewor- 
thy that this method allows all the scattering poles located in 
the domain of interest to be found [6|, [ 16|, [22 1 . Nevertheless, 
these articles made some assumptions on the form of the 
scattering matrix, which usually do not take place for the 
scattering matrix of a diffraction grating. 

In this article we propose new numerical methods for cal- 
culating poles of the scattering matrix with better performance 
(convergence rate, computational complexity, attraction basins 
shape). The paper is organized in six sections. Following the 
introduction, Section [TT] defines the scattering matrix of the 
diffraction grating. In Section [HI] we rigorously derive the 
resonant representation of the scattering matrix with use of 
analytic matrix-valued functions theory. The reader who is 
mainly interested in practical implementation of the proposed 
methods can skip this section. In Sections iTVl HV-Al we review 
the known methods for calculating poles of the scattering 
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matrix. In Section IIV-BI we present new iterative methods 
for calculating poles [Eqs. ( fT9] l, (|20|il, which take account 
of the scattering matrix form in the resonance vicinity. In 
Section IIV-CI a new formulation of the Cauchy-integral-based 
method [Eq. d26l il is proposed. In Section [V] we conduct 
the numerical comparison of the existing methods with the 
proposed ones. 

II. Scattering matrix 

Let us give the definition of the scattering matrix for the case 
of 2D diffraction grating. Assume that the grating is periodic 
along the x-axis, with period d (see Fig. [TJ. Let there be a set 
of plane waves incident on the grating from superstate and 
substrate, with the x-component of the incident light wave- 
vectors given by 



K, 7 



2tt 



m £ Z. 



(1) 



Due to diffraction by the grating, the set of incident plane 
waves is scattered into a set of reflected and transmitted 
diffraction orders. Note that according to Floquet theorem, the 
^-components of the diffraction orders' wave-vectors are also 
defined by Eq. (Q~|i. 

The grating scattering matrix S relates the complex ampli- 
tudes of the incident (vl/ mc ) and scattered (\l> scatt ) waves by 
the formula: 



= * mc , (2) 



where <f lnc = 



h 
h 



R 
T 



R and T are the 



complex amplitude vectors of the reflected and transmitted 
diffraction orders, whereas I\ and I2 are the complex am- 
plitude vectors of the waves incident on the grating from 
superstrate and substrate. Thus, the scattering matrix element 
(S)ij is the complex scattering amplitude of the j-th incident 
wave in the direction of the i-th scattered wave. Note that 
the diffraction of the incident waves by the grating produces 
an infinite number of propagating and evanescent orders, so 
that the matrix S also has infinite dimensions. To conduct 
the numerical simulation, the scattering matrix is truncated. 




Fig. 1. (Color online) Silver diffraction grating on dielectric waveguide-layer 
(grating parameters: period d = 1000 nm, slit depth h gI = 50 nm, slit width 
a = 200 nm, waveguide layer thickness h = 800 nm, layer's permittivity 
e = 5.5). 



Let N — 2K + 1 denote the truncation order, which corre- 
sponds to the incident and scattered waves with the numbers 
m = —K, ■ ■ ■ , K in Eq. (Q~|). Then, the scattering matrix relates 
2N incident and 2N scattered waves. Considering that the 
incident and scattered waves have two polarization states, the 
scattering matrix dimension is AN x AN. The calculation of the 
truncated scattering matrix can be carried out using the Fourier 
modal method [18| also referenced as the Scattering matrix 
method H. In general, the scattering matrix can be calculated 
for multilayer coatings, 2D and 3D periodic structures as well 
as for aperiodic structures (with use of PML or absorbing 
layers). 

For a given geometry and materials of the grating, the 
scattering matrix S is a function of frequency u and the 
x-component of the wave-vector of the incident wave with 
number m = 0: S — S(u),k x ,o)- When describing the 
resonances in the grating transmittance and reflectance spectra, 
the incident wave direction is specified and the scattering 
matrix is treated as a function of frequency u: S = S (oj). 
In this case the structure's resonances are found as the poles 
of the analytic continuation of matrix S (oj) Q. The real part 
of a pole corresponds to the frequency of the incident wave 
that can excite the corresponding eigenmode, while the inverse 
of the imaginary part defines the lifetime of the resonance. 

III. Resonance representation of the scattering 

MATRIX 

In this section we rigorously derive the resonant repre- 
sentation of the scattering matrix © with use of analytic 
matrix-valued functions theory. We believe that this section 
is important, because it gives a rigorous basis for the methods 
proposed in Section [TV] Nevertheless, the reader interested 
mainly in the practical application of the numerical methods 
can skip current section and keeping Eq. © in mind go to 
Section HV] 

Consider the analytic continuation S(oj), oj € C of the 
scattering matrix onto a complex w-plane region D bounded 
by a closed curve T. We assume that the Rayleigh anomalies 
are located far away from the frequency range of interest. 
Then, the analytic continuation S(oj) in the region D will 
be single- valued 0"), ||23l . 

Assume that in the D region there is a simple pole of 
the scattering-matrix analytic continuation at lu = us p . Below, 
only the simple poles of the scattering matrix are considered. 
This assumption is widely used in scattering theory both in 
electrodynamics and quantum mechanics @. In this case, it 
makes sense to define the S(ui) matrix residue: 

Res S(u) = 7~~t S(u) dw, (3) 

w=w p 27T1 J 1 

where the integration contour 7 is chosen in such a way 
as to contain just a single pole u p . Equation ([3]l should be 
understood as an element-wise operation. If the scattering 
matrix has a single pole in the region D, the following relation 
takes place: 

S(u)=A(u) + —^—, (4) 
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where B — Res w=LJ S(ui) and the matrix-valued function 
A(tS) has no poles in the region D, thus being holomorphic 
in this region. In the general case of M poles found in the 
region D, the decomposition in Eq. takes the form: 



M 



B,, 



(m) ' 



(5) 



where B m — Res w w ( m > S(lu). The first term in Eqs. (2) 
and (|5]l describes the non-resonant scattering of light and the 
second — the resonant scattering. 

Let us now discuss the matrix properties of S(lu) in the 
context of scattering (diffraction) theory. Assume that S(lu) 
has a simple pole at lj = lu p , while the inverse matrix S , ~ 1 (w) 
elements have no pole at this frequency. In this case, the 
kernel of the S~ 1 (lo p ) matrix defines non-trivial solutions of 
the homogeneous equation: 



S~ l ijf e 



= 0. 



(6) 



Thus, kci S~ 1 (tjj p ) describes the field distribution in the ab- 
sence of incident waves, i.e. the frequency u p corresponds to 
the quasiguided mode of the grating. 

It can easily be shown that Im B = ker S^ 1 (w p ) and hence, 
rankB = dimkerS' _1 (a;p) (6). As a rule, rankB = 1, 
i.e. only one mode corresponds to the frequency lj p . How- 
ever, at certain grating parameters, the frequencies of several 
different modes may coincide. In this case, rankS > 1. 
Such resonances are referred to as degenerate. The degenerate- 
resonance structures have interesting optical properties which 
can be applied to the design of channel drop filters and all 
pass filters ifTol. ifTTl. 

In the general case, assume that rank B = r. Then, for the 
matrix B E C" XTl , a rank factorization can be written [24]: 



B = LR, 



(7) 



where L G C™ xr , R 6 C rxn , rankL = ranki? = r. In view 
of Eq. ©, Eq. can be represented as 



S = A(lu)+L- 



1 



-R. 



(8) 



Accordingly, the general decomposition in Eq. ((5) takes the 
form: 



S = A(u)+^2 L r , 



1 



rn — 1 



UJ — UJ- 



j^R m = AiujJ+Liluj-n^R, 

(9) 



where L and R are the block matrices defined as L 
[L x L 2 ■■■ L M ], R T = [Rj RJ ••• R T 



and fl p is the diagonal matrix composed of uj p m \ m = 
1, M, with the frequency Wp™' recurring as many times as 
is the rank of the corresponding matrix B m . Equations (|S), (O 
define resonance representations of the scattering matrix which 
will be used in further developments of numerical methods for 
computing the scattering matrix poles. 



M 



IV. Computing the scattering matrix poles 

As demonstrated above, the poles of the scattering ma- 
trix define the grating eigenmodes. Consider a problem of 
computing the scattering matrix S(ui) poles located in the 
region D. The simplest approach assumes that the poles of 
the matrix-valued function S^u;) can be found as the poles 
of its determinant, thus by solving numerically the following 
equation: 

l/detS(w) = 0. (10) 

This approach allows the poles of a small-dimension matrix 
S(uj) to be found. For a large-dimension matrix S(uj), the 
calculation of the determinant becomes numerically unstable. 
Showing a higher stability is the problem of solving the 
equation [16|: 

l/maxeigS(u;) = 0, (11) 

where maxeigS^w) is the maximal-modulus eigenvalue of 
the matrix S(uj). Equations (TlOb and (fTTT i can be solved 
with iterative methods for deriving the root of a nonlinear 
equation, for instance, Newton's method or, more general, 
Householder's method ll25l . For Eq. (fill , this method is given 
by the following iterative procedure: 



d" 



0J n+ X = 0J n + p 



— maxeig S(co) 



g^maxeigSV) 



(12) 



where cj„ is the initial pole guess. At p = 1, Eq. (fT2l) 
corresponds to Newton's method; at p = 2 — to Halley's 
method. When solving Eqs. ( TTOb and (fTTT l. the method of 
Eq. (fT2l disregards the form of the scattering matrix in Eq. (O 
replacing it with a scalar value [determinant or maximal- 
modulus eigenvalue of S(ui)]. This is the reason of rather 
slow convergence of the above mentioned methods. In the 
following subsections we consider the iterative methods for 
computing the scattering matrix poles with the use of matrix 
decompositions. 



A. Computing the poles through linearization of the scattering 
matrix inverse 

Let us analyze the iterative method for computing the 
scattering matrix poles proposed in Refs. Q, ||8], lfP7l . [21 1. 
Let uj n be an initial pole guess. We decompose the matrix 
S , ~ 1 (cj) into a Taylor series up to the first term: 



S- 1 (u) = S- 1 (u n ) + 



dS- 



duj 



(13) 



Assume that u; p is the scattering matrix pole, then there exists 
a vector \|/ scatt that represents a non-trivial solution of the 
system ©. Multiplying Eq. ( fl3l l by \J/ scatt on the right at 
oj = lu p yields: 



(uj n - u p ) 



dS- 



(14) 



Equation (TT~4T > defines a generalized eigenvalue problem. Solv- 
ing this problem yields a set of eigenvalues = ui n — lu p . 
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Choosing out of these a minimal-modulus eigenvalue, we 
obtain the following iterative method: 



u n+ i = Lu n - mineig S 1 (u n ), 



dS- 



(15) 



where eig(F, G) denotes a vector composed of the eigenvalues 
for the generalized eigenvalue problem, FX — A GX. Choos- 
ing in Eq. ( fT5l l a minimal-modulus eigenvalue means that the 
subsequent guess io n+ i chosen will be closest to the initial 
guess oj n . 

Note that the method of Eq. ([T5T > has a number of disad- 
vantages, the major being inability to conduct the numerically 
stable computation of the scattering matrix inverse for a large 
number of diffraction orders [21 ]. Another disadvantage is that 
if a pole of the matrix is located in the vicinity of the 

pole of the scattering matrix S(uS), the Taylor series (IT3b will 
diverge (see Appendix IBl. 

B. Pole calculation based on a resonance approximation 

Let us consider a new iterative technique for scattering 
matrix pole calculation that takes account of the scattering 
matrix form in the resonance vicinity. Assume that an initial 
pole guess is cj = uj n . Put down the resonance approxi- 
mation © neglecting the frequency-dependence of the non- 
resonance term: 



S(u) = A + L{ul - ilp^R. 
The first and second derivatives of S(cj) are 

s'(w) = -L{ui- n p y 2 R, 

S"(u>) = 2L(uI - il p )- 3 R. 



(16) 



(17) 



Assume that ranki = rank/? = ^ TO rankL m = 
2~^ m ranki? m . This assumption implies that the columns of 
the matrix L are linearly independent or, which is the same, 
the kernels of the matrices S -1 ^^), m = 1,...,M are 
linearly independent. The latter means that the scattered field 
distributions for different modes are linearly independent. 
This assumption normally holds when the scattering matrix 
dimension is much larger than the number of modes (dim S > 
rankL). 

Putting lu = uj n in Eq. ( fTTT i yields a system of two matrix 
equations with respect to an unknown diagonal matrix f2 p . 
A method for solving equations of this kind is considered 
in Appendix [A] Following that method [see Eq. OTTll. the 
diagonal matrix ft p is given by 



r! p = w „/+2diageig(C/ r t 5'K)T4S; 1 ), 



(18) 



where diageig F is the diagonal matrix composed of the 
eigenvalues of the matrix F and the matrices U r , E r , V r are 
derived from the compact singular value decomposition of the 
matrix S"(uj n ) = U r 'E r VJ. From Eq. (fT8l . we can obtain the 
following iterative method: 

w„+i = uj n + 2 mineig (l# S' '(w^X^ 1 ) . (19) 

In a similar way to Eq. $15[ , taking the eigenvalue with the 
minimal modulus in Eq. ( fT9l denotes that the subsequent pole 



guess chosen u) n +i will be the pole closest to the initial guess 

The iterative method in Eq. ( TT9i l assumes that there are 
several poles in the initial guess vicinity. In practice, one 
can assume that in the vicinity of Lu n there is a single pole 
corresponding to the non-degenerate resonance. In this case, 
in view of Eq. d34l i, the iterative method takes a simple form: 

maxeig S'(uj n ) 



oj n+ i — u) n + 2- 



(20) 



maxeig S"(u n ) 

As distinct from ([T5l l. the iterative methods in Eqs. JT9l 
and (f20b are based on the resonance approximation ([TBI , rather 
than the linearization of the matrix S ,_1 (o;). The results of 
the numerical investigation in Section [V] demonstrate that the 
iterative method $1% and (l20l show a better convergence 
when compared with Eq. (fT5T >. Besides, the benefit of the 
proposed approach is that it remains valid for a large number 
of diffraction orders (large-dimension matrix S) and when the 
poles of the scattering matrix and of its inverse are close to 
each other. 

Note that the method of Eq. (fTst has the meaning of 
Newton's method for the matrix-valued functions, whereas 
the methods of Eqs. ( fT9b and d20i > can be treated as a matrix 
extension of Halley's method for solving equations of the form 
1/ f(x) — 0. In general, it is possible to write down a matrix 
analog of Householder's method |25| with the aid of the p-th 
and (p— l)-th derivatives of the scattering matrix: Sw(w„), 
S(P- 1 )(w n ). In this case, the iterative method of Eq. ( fT9l will 
represent a particular case (at p — 2) of the following iterative 
method: 



w n+ i =oj n +p mm eig 



Vr'Sr 



(21) 



where the matrices U r ,'S r ,V r are derived from the com- 
pact singular value decomposition of the matrix S^ pS> (uj n ) = 
U r Y, r V^ . The analog of the method in Eq. (l20l is written in 
a similar way: 

maxeig SO" 1 )^) 



P- 



(22) 



maxeig S^(uj n ) 

The above relation resembles Householder's method of 
Eq. (ITZb for solving Eq. (fTTT i. except that the operations 
of differentiation and calculation of the maximal-modulus 
eigenvalue are swapped. At p = 1, the methods in Eq. ( |2"TT > 
and d22i > represent analogs of Newton's method, but, as distinct 
from Eq. $15[ , the said method relies upon the calculation of 
the scattering matrix S(uj) and its derivative, rather than of 
matrix 

C. Pole calculations based on the Cauchy integral 

The efficiency of the above-discussed methods in 
Eqs. (|T0l>, CO}, CSJ, (HU, and (|20j essentially depends on 
the initial pole guess lu = cj n . Besides, the said methods 
allow one to find just a single pole in the initial guess' 
vicinity. Note that with methods of Eqs. (fT5l ). ( fT9l ) and ( T2TI ) 
it is possible to use several eigenvalues when calculating the 
subsequent approximations. However, such an approach does 
not guarantee that all the poles of the scattering matrix will 
be found in the D region of interest. 
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Below, we discuss an approach based on the Cauchy integral 
that enables all the scattering matrix poles to be found in the 
region of interest. Such an approach is widely utilized when 
calculating modes of multilayered slab waveguides [2|, [3| and 
photonic-crystal waveguides [6|, [16|. 

Following Refs. (6), fl6l . ||22| . calculate two contour inte- 
grals with respect to the matrix-valued functions: 



C 2 = i> uS(uj) dui, 

Z7T1 



(23) 



where T is the boundary of the region D. Making use of 
Eq. dD, we easily find that 



C\ = B rn ; 



C 2 



(m) 



Bn 



(24) 



where B m = Res = (m) S(uj). Note that in Ref. it 
was demonstrated with a number of examples that all the 
poles found within the contour T could be derived using the 
integrals (l23l . It should be noted that in Ref. |6| the matrices 
B m are assumed to be symmetrical and the columns of the 
matrix L in Eq. © — orthogonal. The said assumptions 
are usually not valid for the diffraction grating scattering 
matrices (8]. 

Let us discuss the method for solving Eq. (l24l > under more 
general assumptions. To these ends, Eq. (f24-b is written using 
the notations of Eq. ©: 



C\ — LR', 
C 2 = LQ.pR. 



(25) 



As earlier, the columns of the matrix L are assumed to be 
linearly independent. This condition is more general than the 
orthogonality condition employed in Ref. [6|. The system of 
equations (l25l > will be solved with respect to the unknown 
diagonal matrix 51 p using a method discussed in Appendix lAl 
Following this method the Sl p matrix is given by 



ftp^diageig^C^E- 1 ) 



(26) 



where U r , S rj V r are derived from the compact singular 
value decomposition of the matrix C\ = f/ r E r V^J. 

Thus, calculating the integrals in Eq. ( l23l enables all the 
scattering matrix poles located inside the integration contour 
to be found. The pole calculation accuracy depends on the 
accuracy of the numerical calculation of the integrals (123V 
To calculate the poles with high accuracy, a large number of 
the scattering matrix calculations needs to be conducted [6|, 
iTTBI . Because of this, it is relevant to use the Cauchy-integral- 
based method only when calculating the initial pole guess, 
with the subsequent highly accurate values to be derived using 
the above-considered iterative techniques. 

V. Numerical examples 

In this section we study the performances of the above- 
discussed methods for computing the poles of the diffraction 



grating scattering matrix. The diffraction grating scattering 
matrix is calculated using the Fourier modal method [18], [ 19|, 
[26|. Note that the condition of analyticity of the scattering 
matrix requires proper implementation of the Fourier modal 
method in the case of complex frequencies 101, ifTTl . 

Let us consider the calculation of the eigenmodes of the 
diffraction grating composed of a periodic array of slits in a 
silver film coated on a dielectric waveguide layer (Fig. [T). We 
are interested in the modes that can be excited by the normally 
incident plane wave (fc x> o = 0). The grating parameters are 
given in the caption to Fig. Q] The permittivity of silver is 
described using the Lorentz-Drude model ||271 . with the an- 
alytical approximation for £A g (w) being treated as a function 
of complex frequency. Note that the analyticity of the function 
e(ui) is required to provide the analyticity of the scattering ma- 
trix S(lj). If this condition is violated, the convergence rate of 
the iterative techniques will be essentially deteriorated, while 
the Cauchy-integral-based method will become inapplicable. 
Below, the grating mode frequency ui p is assumed to have been 
calculated correctly if maxsvd S(u p ) > 10 10 , where maxsvd 
denotes maximal singular value. This condition is also used 
as a termination criterion for iterative methods. The numerical 
calculations show that for the grating under study the above 
condition corresponds to the pole determination accuracy of 
AAp| < nm. Such an accuracy is sufficient for practical 
uses, including the calculation of the mode field distribution. 

A. Cauchy-integral-based method 

Let us determine the entire set of modes supported by the 
grating in the frequency range to 6 [1.0 x 10 15 s ; 1.8 x 
10 15 s -1 ]. Note that only high quality modes are considered 
(|Imcjp| < 5 x 10 12 s _1 ). We used the value of N = 21 as 
the truncation order for subsequent calculations. 

For estimating the mode frequencies, we calculated the 
contour integrals of Eq. (1231 using the trapezoidal method 
with 500 discretization points. Ten approximate pole values 
were derived from Eq. (|26V The said values are shown in col- 
umn 2 of Table U The accurate values of the pole frequencies 
calculated by the iterative method of Eq. (f20]l are shown in 
column 3 of Table |T] The initial guess error is presented in 
column 4. Table U suggests that with the contour integration 
method, the poles can be calculated with a sufficiently high 
accuracy. 

Note, for comparison, that if the pole approximations were 
obtained by calculating the scattering matrix on a 500 x 20 grid 
(i.e. 10000 calculations of the scattering matrix), we would be 
able to find only nine poles. The reason is that the poles ujq, uj-j 
are located very close to each other. Thus, the use of the 
Cauchy-integral-based method is appropriate in this case. 

B. Comparison of the convergence rates of the iterative meth- 
ods 

Consider the convergence of the iterative methods of 
Eqs. ( fTBI l, d20l l. and Newton's method when solving Eq. dTQb - 
First, analyze the situation when fairly good initial guesses of 
the poles are known. In this case, the pole estimates shown 
in column 2 of Table Q] can be refined using the methods 
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Poles derived from Eq. ( |26t and a more accurate iterative pole calculation. 



niter 





0>Cauchy 


^-"accurate 


Inaccurate 0?Cauchy | 


. 

Eq. (20) 


. . 

Eq. (15) 


Newton 
for 

Eq. (10( 




1.14246 X 10 Ib 


- 3.35055 x lO^i 


1.14247 X 10 lb 


- 3.34402 x 10 I2 i 


7.86650xl0 M 


1 


2 


3 


W2 


1.15554 x 10 15 


- 6.85407 x 10 n i 


1.15551 x 10 15 


- 6.86647 x 10 n i 


3.40401 xlO 10 


2 


2 


3 


0)3 


1.40431 x 10 15 


- 1.87394 x 10 12 i 


1.40428 x 10 15 


- 1.87643 x 10 12 i 


3.45275X10 10 


2 


2 


3 


U?4 


1.48168 x 10 15 


- 3.86736 x 10 12 i 


1.48168 x 10 15 


- 3.86743 x 10 12 i 


2.35895 xlO 9 


1 


2 


3 


<*>5 


1.50470 x 10 15 


- 1.11064 x 10 12 i 


1.50469 x 10 15 


- 1.11116 x 10 12 i 


1.52001X10 10 


2 


2 


3 


UJ6 


1.66425 x 10 15 


- 7.11868 x 10 n i 


1.66425 x 10 15 


- 7.12688 x 10 n i 


3.44210 xlO 9 


1 


2 


3 


0>7 


1.66551 x 10 15 


- 1.05160 x 10 n i 


1.66545 x 10 15 


- 1.06986 x 10 n i 


5.64199X10 10 


2 


2 


4 


0>8 


1.73171 x 10 15 


- 3.70973 x 10 12 i 


1.73171 x 10 15 


- 3.70850 x 10 12 i 


1.23369 xlO 9 


1 


2 


3 


U>9 


1.73768 x 10 15 


- 2.43030 x 10 12 i 


1.73768 x 10 15 


- 2.43055 x 10 12 i 


2. 84178 xlO 9 


1 


2 


3 


0J10 


1.77660 x 10 15 


- 2.79575 x 10 12 i 


1.77659 x 10 15 


- 2.79747 x 10 12 i 


6. 09846 xlO 9 


1 


2 


3 



under study. Columns 5, 6, 7 of Table U give the number of 
iterations when calculating the poles by the methods (l20l i. ( [T5| >. 
and Newton's method when solving Eq. ( TTOb . From Table U 
it follows that the convergence rate is lowest for Newton's 
method when solving Eq. ( TTOb . Meanwhile, the proposed 
method of Eq. ( f20b and the method of Eq. dT3T > require about 
the same number of iterations. 

Let us investigate in more detail the convergence rate when 
the initial pole guesses are not known. To do so, the following 
numerical experiment was conducted. For a random point from 
the complex frequency range under study, the iterative method 
was initiated. Tables [TTJ [TTTl [IV] give the average percentage of 
the initial guesses having converged to the pole (n conv ), the 
average number of iterations (niter), and the average number 
of the scattering matrix calculations (n s _ ca i c ). 

Note that different methods require a different number of 
the scattering matrix calculations per iteration. For instance, 
Newton's method and the method of Eq. ( fT3l ) calculate the first 
derivative of the scattering matrix, requiring two scattering ma- 
trix calculations per iteration, whereas Halley's method and the 
methods of Eqs. ( TT9i l. d20l i proposed here calculate the second 
derivative, requiring three scattering matrix calculations per 
iteration. Therefore, for the comparison purposes, what should 
be taken into consideration in the first place is the magnitude 
of n s - C aic that defines the total number of the scattering matrix 
calculations. At the same time, in each iteration, the scattering 
matrices can be calculated independently, or in parallel. Thus, 
for parallel implementation of the iterative methods, it is the 
magnitude of niter which turns out to be of greater importance. 

Table llll was calculated for the truncation order N = 21. It 
can be seen from Table [II] that among the "scalar" iterative 
methods for solving Eqs. ( TTOb and ( TTT1 >. Halley's method 
requires the smallest number of iterations. This is due to 
the fact that Halley's method has cubic convergence, whereas 
Newton's method converges with order two 11251 . Note also 
that Newton's method shows better convergence for Eq. (ITOl 
than for Eq. (ITTb . 

In terms of the number of the scattering matrix calculations, 
the matrix methods of Eqs. (I151 l. < fT9b . < f20b possess a higher 
rate of convergence when compared with the scalar methods. It 
is seen from Table HIl that the methods ([T9l i and (l20l > proposed 
here are able to find the poles within a smaller number of 
iterations among all the methods discussed in this paper. In 



TABLE II 

Convergence of the iterative methods for the truncation 
order TV = 21 (84 x 84 scattering matrix). 



Method 


f^conv t' 


"iter 


^s — calc 




Newton for Eq. ill 01" 


98.5 


9.85 


19.7 


57.2 


Newton for Eq. ill U" 


94.9 


7.08 


14.16 


42.4 


Halley for Eq. QTJ" 


99.9 


4.20 


12.60 


54.5 


Eq. (L5} b 


99.7 


5.48 


10.96 


63.9 


Eq. (I9) b 


100 


2.78 


8.35 


96.1 


Eq. HD b 


100 


3.40 


10.2 


64.6 



"Scalar method 
''Matrix method 

TABLE III 

Convergence of the iterative methods for the truncation 
order tv = 61 (244 x 244 scattering matrix). 



Method 


^conv 


n itcr 


^•s — calc 


'^ncar 


Newton for Eq. (10) 


82.5 


11.79 


23.58 


37.7 


Newton for Eq. ill 11 


87.2 


7.1 


14.2 


24.5 


Halley for Eq. QT} 


99.7 


4.54 


13.62 


40.3 


Eq. {B) 


78 


5.95 


11.9 


48.6 


Eq. 03 


100 


2.88 


8.64 


87.6 


Eq. |20) 


100 


3.72 


11.16 


43.75 



TABLE IV 

Convergence of the iterative methods for the truncation 
order TV = 201 (804 x 804 scattering matrix). 



Method 


flconv 


,% n itcr 


^s — calc 




Newton for Eq. (101 











Newton for Eq. (Ill 


91 


6.83 


13.66 


26.5 


Halley for Eq. 0TJ 


99 


4.85 


14.55 


40 


Eq. (B) 











Eq. (19) 


100 


3.86 


11.58 


75.5 


Eq. {20) 


100 


3.94 


11.82 


45 



particular, the average number of iterations for the method (fT~9b 
is nearly twice as small as that for the method ([T5V At 
the same time, the average number of iterations required to 
calculate the scattering matrix by the method ([19} is just 1.31 
times smaller than for the method of Eq. (Q3). The comparison 
of the methods (fT9b and ( f20b shows that the latter approach 
based on the assumption of a single resonance in the vicinity 
of the initial pole guess requires on the average a 1.22 times 
larger number of iterations than a more general method ([19} . 

Let us analyze how the methods under discussion operate at 
a larger truncation order. Table lUsuggests that at TV = 61 the 
convergence of the method dT31 > and Newton's method when 
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solving Eq. (fTOb is deteriorated: with the pole finding rate 
being 78% for the former method and 83% for the latter. When 
N = 201 (Table |IV]), both the method (Q~5]) and Newton's 
method for solving Eq. ( fTOb become inoperative. When used 
to solve Eq. ( TTTb . Newton's and Halley's methods are converg- 
ing, although requiring by about 30% more iterations when 
compared with the methods f!9[ and (l20l >. The same is true 
regarding for the number of the scattering matrix calculations. 
It is noteworthy that the methods ( fT9b and (l20b converge in 
one hundred percent of the cases. 

In addition to the rate of convergence, an important char- 
acteristic of iterative methods is the shape of the basins of 
attraction. With regard to the calculation of the poles, basin 
of attraction of the pole lj p is the set of points from which 
the iterative method converges to ui p . Ideally, one would like 
the iterative method to converge to the nearest pole. In this 
case, if we vary some parameter of the grating, we will be 
able to calculate the corresponding variation of the mode 
eigenfrequency. Hence, we can easily calculate the functional 
dependence of the eigenfrequency on some parameter (for 
example we can calculate the dispersion curves by varying 
the wavenumber k x ). Unfortunately, the Newton's method 
(as well as the Halley's method) is characterized by the fractal- 
shaped basins of attraction |28l . 

Fig. |2] shows the basins of attraction for Recjo G [1.65 x 
10 15 s _1 ;1.75 x 10 15 a -1 ], Imw G [-5 x 10 12 s^ 1 ^ x 
10 12 s" 1 ] region. We used the value of N = 21 as the 
truncation order. According to Fig. |2] the basins of attraction 
for the proposed method ( [T9| > have rather regular shape while 
other methods give pronounced fractal shape. 

As the numerical characteristic of the shape of the basins of 
attraction consider the percentage of the starting points from 
which the method converges to the nearest pole. This value 
(rineai ) presented in the last column of Tables IIThIVI has to be 
close as possible to 100%. As follows from Tables HIhIVI the 
proposed method ( [T9| > in the most cases (70-90%) converges 
to the nearest pole. At the same time for the other methods 
the value of n ncar ranges from 20 to 60%. 

Tables HH4IVI demonstrate that the performance of the con- 
sidered methods decreases with increasing the dimension of 
the scattering matrix. It must be noted that in all above- 
mentioned methods except method ( fT5l ) it is possible to use 
only a part of the large-dimension scattering matrix. For exam- 
ple, we can use the central part of size 84 x 84 from the large- 
dimension scattering matrix computed at the truncation order 
N = 201. In this case, it is expected that the performances 
of the matrix methods of Eqs. (1% , (|20| | as well as Newton's 
(Halley's) method when solving Eqs. ( [Tol l. (fTTT i will be close 
to the values given in Table HU 

VI. Conclusion 

Summing up, a new iterative technique for calculating the 
scattering matrix poles has been proposed in this work. Taking 
account of the scattering matrix form in the pole vicinity, the 
method relies upon solving the matrix equations through the 
use of a matrix decompositions. Using the same mathematical 
approach, a Cauchy-integral-based method that enables all 



the poles in a designed region to be calculated has been 
described. The calculation results for the modes supported 
by metal-dielectric diffraction grating have shown that the 
iterative method proposed has the high rate of convergence and 
is numerically stable for large-dimension scattering matrices. 
The proposed method converges to 30-60 % faster compared 
to known methods. An important advantage of the proposed 
method is that it converges (in the most cases) to the nearest 
pole on the complex plane. 

Note that the methods discussed in the present article treat 
the scattering matrix as a function of complex frequency u. 
However, all the methods remain valid if the scattering matrix 
is treated as a function of the incident light wave number 
k x ,o 13 or as a function of the z-component of the wave- 
vector of a diffraction order (k z ,„). In the case of complex 
ui, the method can be used to explain the resonances in the 
transmittance (reflectance) spectrum of the structure. The use 
of complex k x $ makes the proposed methods an effective tool 
to calculate the guided or leaky modes of the grating waveg- 
uides. The complex k Zjn approach is employed to describe the 
modes in the vicinity of Rayleigh frequencies of the diffraction 
grating fll, ES- 

Moreover, the considered methods are not limited to the 
optical scattering but also can be applied to the problems of 
acoustical and quantum mechanical scattering. 

Appendix A 
Solving a system of matrix equations 
A = LXR, B = LYR 

Given matrices A,B€ C nx ™ that take the form: 



A = LXR; 
B = LYR, 



(27) 



rxn are some unknown matrices, 
are unknown diagonal matrices, 



where L e C nxr , R G 
X G C rxr , and Y G C rxr 
with the value of r being also unknown. It is required to find 
the relation between the matrices X and Y. Note that systems 
of this type are encountered when finding the poles of Linear 
Time Invariant (LTI) systems 11291, QUI. 

Consider a method for solving the system (l27t which is 
based on the singular value decomposition 11311 . First, we 
find the rank r of the matrix A. To do so, we calculate its 
singular value decomposition: A = UT,V^ (here V"f denotes 
the conjugate transpose of V). The rank r is defined by the 
amount of nonzero singular numbers. In practice, it may occur 
that the A matrix is noised. In this case, r is chosen to be equal 
to the amount of singular numbers larger than a threshold 
value. Knowing the rank of A, the compact singular value 
decomposition is given by 



A=U r T, r Vj, 



(28) 



where S r G R rxr is a diagonal matrix composed of r 
largest singular numbers of the A matrix; U r ,V r G C nxr 
are the matrices of the corresponding left and right singular 
vectors. If the matrix is noised, the equality in Eq. (l28l 
becomes approximate. In this case Eq. d28l represents the 
rank-r approximation of matrix A. 
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(29) 



Note that both U and V are unitary, hence U r = VJV r — 
I. Multiplying the equations in d27b by C7J on the left and by 
V r on the right, we obtain: 

(U}L)X(RV r ) = {W r U r )Y, r {V}V r ) = S r 
{W r L)Y{RV r ) = W r BV r . 

Note that rank£ r = rank A = rankC/ r = rankV,. = r, 
hence both UJL and RV r are square and invertible matrices 
from C rxr . Multiplying the second equation in d29b by E^T 1 
on the right, we get: 

{U^YX^iUlL)- 1 = UlBVrY,-; 1 . (30) 

Equation (f30b can be treated as the eigendecomposition of the 
matrix on the right-hand side of the equation. Thus, the YX^ 1 
matrix can be explicitly given by 

YX- 1 = diageig([/ r t W r E- 1 ). (31) 

Consider an important particular case when the matrices A 
and B are of unit rank, i.e. r = 1. In this case, the vector 
Ui is the eigenvector of the matrix A and the vector V\ is 
the eigenvector of the matrix A*. Hence, the A matrix can be 
represented as 

maxeig A^ 



A = U-l- 



(32) 



where max eig A is the only non-zero eigenvalue of the matrix 
A. The vectors U\, Vi are also the eigenvectors of the matrices 
B and B\ respectively. Thus, 

maxeig B 



B = Ux- 



Substituting (O and (O into (EB yields: 

maxeig B 

y/x 



(33) 



(34) 



maxeig A 

Note that because in practice the matrices A and B can be 
noised, the maximal-modulus eigenvalues should be used in 
Eq. (01. 

Appendix B 

An example of closely located poles of the 

MATRICES S(uj) AND S* -1 ^) 

In this appendix, we discuss a specific diffraction grating for 
which the iterative method in Eq. ([l5T l fails to converge already 
at a small value of the truncation order. Let us consider a silver 
diffraction grating (period d = 1000 nm, height h gI = 10 nm, 
slit width a = 990 nm) on a silver substrate. 

Analyze how the method of Eq. ( fl3] l performs for the 
truncation order N = 21. We are interested in the pole 
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with the frequency lo p w 1.859 x 10 15 - 1.613 x 10 12 s _1 . 
The peculiarity of this structure is that the scattering matrix 
determinant turns to zero in the vicinity of ui p , at lu rs 



1.863 x 10 15 - 1.465 x 10 12 s" 



Let the initial guess be uj — 1.86 x 10 5 s . With this ini- 
tial guess, the proposed iterative method (|20T > converges within 
two iterations; whereas the iterative method (TTSt converges 
within 5 iterations. However, for a bit worse initial guess of 
ujo = 1.85xl0 15 s _1 , while the proposed iterative method (l20l 
is still converging, the method ([T5l l fails to converge. The 
reason why the method (TT3T > fails to converge is that 5~ 1 (cli) 
has a pole and, consequently, the convergence radius of its 
Taylor series in Eq. JTST l is small. 
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